********* Correlating geographic subsidy, CZ F.E., price index, and log median household income 

	use "${intermediate_data}/twfe/geographic_variation_cz90level.dta", clear
	
	keep if !missing(fe_cz_physicians)
	
	foreach var of varlist GAF a0_laspeyres med_hhinc2016  {
		gen ln_`var' = ln(`var')
	}
	gen diff_GAF_laspeyres = ln_GAF - ln_a0_laspeyres

	eststo clear

	eststo: reg diff_GAF_laspeyres ln_med_hhinc2016, robust
	g elast_sub_inc = e(b)[1,1]
	g n_ur_sub_inc = e(N)

	eststo: reg ln_GAF ln_med_hhinc2016, robust
	g elast_gaf_inc = e(b)[1,1]
	g n_ur_gaf_inc = e(N)

	eststo: reg ln_a0_laspeyres ln_med_hhinc2016, robust
	g elast_prices_inc = e(b)[1,1]
	g n_ur_prices_inc = e(N)
	
	eststo: reg fe_cz_physicians ln_med_hhinc2016, robust
	g elast_czfephysicians = e(b)[1,1]
	g n_ur_czfephysicians = e(N)

	*Output elasticity estimates
	esttab using "${temp}/govtpolicy01-cz_gaf_elasticity.csv", replace cells(b(fmt(%09.4g)) se(fmt(%09.4g))) ///
	constant starlevels( * 0.10 ** 0.05 *** 0.010) stats(N, labels("N_UR")) ///
	mtitles("$\text{Log}\left(\dfrac{\text{GAF}}{\text{Price Index}}\right)$" ///
		"Log GAF" ///
		"\shortstacking{Diamond and Moretti (2021) \\ Price Index}" ///
		"\shortstacking{CZ Fixed Effects \\ Physicians}") /// 
	coeflab(_cons "Constant" ln_med_hhinc2016 "Log Median Household Income")
	
	* Format for export 
	import delimited "${temp}/govtpolicy01-cz_gaf_elasticity.csv", clear asdoub
	insobs 1, after(7)
	
	ds v1, not
	foreach var of varl `r(varlist)' {
		loc N_UR =  subinstr(subinstr(`"`=`var'[9]'"',"=","",.),`"""',"",.)
		drbcount `N_UR' 
		replace `var' = `"="`r(rounded)'""' if `var' == ""
	}
	replace v1 = `"="N""' if v1 == ""
	
	export delimited using "${mypath}/intermediate_csv/govtpolicy01-cz_gaf_elasticity.csv", replace dataf delim(tab) 

********* Visualizations

	use "${intermediate_data}/twfe/geographic_variation_cz90level.dta", clear
	
	foreach var of varlist GAF a0_laspeyres med_hhinc2016  {
		gen ln_`var' = ln(`var')
	}
	gen diff_GAF_laspeyres = ln_GAF - ln_a0_laspeyres

 	*** Diff_GAF_laspeyres vs. median hh-income

	preserve
		binscatter diff_GAF_laspeyres ln_med_hhinc2016, reportreg genxq(bin)

		reg diff_GAF_laspeyres ln_med_hhinc2016, robust
		g hhinc_cons = _b[_cons]
		g hhinc_slope = _b[ln_med_hhinc2016]
		g hhinc_slope_se = _se[ln_med_hhinc2016]
		g hhinc_r2 = e(r2)
		sum ln_med_hhinc2016 if e(sample)
		gen xmean = `r(mean)'

		gcollapse (mean) diff_GAF_laspeyres ln_med_hhinc2016 hhinc_cons hhinc_slope hhinc_slope_se hhinc_r2 xmean (sum) N_UR = N_UR_physicians, by(bin)
		drop if bin == . 
	
		drbest_docinc fe_cz_physicians diff_GAF_laspeyres fe_cz_diff_cons fe_cz_diff_slope xmean fe_cz_diff_r2
		order N_UR*, last
		export delimited using "$mypath/intermediate_csv/govtpolicy01-binscatter_gaf_diff_med_hhinc2016.csv", replace dataf delim(tab)  
	restore 

	*** Binscatter FEs vs. [ln(GAF) - ln(CZ Prices)]
	
	keep if !missing(fe_cz_physicians)
	
	binscatter fe_cz_physicians diff_GAF_laspeyres, reportreg genxq(bin)

		reg fe_cz_physicians diff_GAF_laspeyres, robust
		g fe_cz_diff_cons = _b[_cons]
		g fe_cz_diff_slope = _b[diff_GAF_laspeyres]
		g fe_cz_diff_slope_se = _se[diff_GAF_laspeyres]
		g fe_cz_diff_r2 = e(r2)
		sum diff_GAF_laspeyres if e(sample)
		gen xmean=`r(mean)'	
		
		gcollapse (mean) fe_cz_physicians diff_GAF_laspeyres fe_cz_diff_cons fe_cz_diff_slope fe_cz_diff_slope_se xmean fe_cz_diff_r2 (sum) N_UR = N_UR_physicians, by(bin)
	
		drop if bin == .
	
		drbest_docinc fe_cz_physicians diff_GAF_laspeyres fe_cz_diff_cons fe_cz_diff_slope xmean fe_cz_diff_r2
		order N_UR*, last
		export delimited using "$mypath/intermediate_csv/govtpolicy01-binscatter_fe_cz_gaf_diff.csv", replace dataf delim(tab)  
